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I. INTRODUCTION 

The problem of the merger of binary black hole sys- 
tems now seems to be well under the control of numerical 
relativity. More precisely, the development, due to Ein- 
stein's vacuum evolution equations, of an initial data set 
containing two apparent horizons into a quiescent state 
containing only one apparent horizon, has now been car- 
ried out numerous times by various research groups, with 
somewhat different numerical treatments and mathemat- 
ical formalisms [H~Q]. Numerical relativity is now a tool 
for studying the physics of strong gravitational hclds. 

When applying this tool, one is immediately faced with 
a fundamental irony of numerical relativity: a numerical 
code is incapable of dealing with abstract tensors, and 
must instead compute their components in a particular 
vector basis. The fundamental physics of general relativ- 
ity, however, is basis independent. One must be careful 
to ensure that any physical claims are independent (to 
whatever extent is possible) of the coordinate system and 
vector basis in which they are demonstrated. 

One reasonably well-developed example is the compu- 
tation of spin angular momentum in binary black hole 
simulations. Numerous investigations have been made of 
the physics of spinning black hole mergers, presenting in 
some detail effects such as a hang-up of the merger, allow- 
ing angular momentum to be radiated so that the final 
remnant has sub-extremal spin [f| ; spin flips @ , in which 
the dynamics of the merger cause the spin direction of the 
merged black hole to be dominated by the direction of 
orbital angular momentum, rather than the spins of the 
progenitor black holes; and perhaps of most astrophysical 
interest, the kick applied to a merged black hole system, 
balancing the linear momentum given off in gravitational 
radiation during nonsymmetric mergers [7H2C| . A certain 
amount of investigation has also gone into the study of 
black holes of nearly-extremal spin in binary configura- 
tions, an avenue that could probe the limits of cosmic 
censorship [2l|, [22j . Because such physical effects must 
be parametrized according to the spin angular momenta 



of the dynamical black holes, methods must be devised to 
define and compute such a quantity. The most common 
approach begins with a formula that appears both in the 
quasilocal formalism of Brown and York |23| and in the 
isolated and dynamical horizon formalisms [24 l25j . This 
formula gives angular momentum within a two-surface 
(normally taken to be an apparent horizon of spherical 
topology) as a functional of a vector field tangent to that 
surface. This vector field is interpreted as a generalized 
rotation generator, and it is through this that the vec- 
torial nature of angular momentum in Newtonian me- 
chanics is generalized. In order to apply this formula, a 
rule must be given for choosing such a generalized rota- 
tion generator on a dynamical black hole. Methods have 
recently been presented to fix these vector fields as "ap- 
proximate Killing vectors" in a precise sens e |22| . l2q - [29j . 

The method presented in Refs. [22|, [28|, l29l] actually 
provides much more information than just the general- 
ized rotation generators. The method starts with the 
expression of the vector field in terms of a scalar poten- 
tial: 



AB 



(1) 



where uppercase latin letters index the tangent bundle to 
the two-dimensional surface, V is the covariant derivative 
on this tangent bundle, inherited from that on spacetime, 
and cab is the Levi-Civita tensor on the surface. The 
vector <fi is said to be an approximate Killing vector if 
it is of this form and if the function z satisfies a certain 
generalized eigenvalue problem on the surface. On a met- 
ric sphere 1 , the operator in this problem reduces to the 
conventional spherical Laplacian, so these functions can 
be interpreted as spherical harmonics of the two-surface. 
In this special case, the three I = 1 harmonics provide 
the three standard rotation generators. 

The appearance of generalized spherical harmonics in 
this formalism raises the possibility that one could nat- 
urally define more than just the spin angular momen- 



1 Throughout this paper, by "metric sphere" we mean a sphere in 
the metric sense: a closed 2-surface of constant positive intrinsic 
'Electronic address: owcn@astro.corncll.edu curvature, sometimes also referred to as a "round sphere." 
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turn (which is often physically understood as the current 
dipole moment of the source). Perhaps with the help 
of the remaining eigenf unctions, we could define higher 
multipole moments. 

The idea of quasilocal source multipoles in general rel- 
ativity is not new. In Ref. (3C| , a complete formalism was 
presented for application on axisymmetric isolated hori- 
zons. This formalism involves numbers /„ and L n , where 
n is a nonnegative integer index. Ashtekar et al. not 
only provided definitions for these multipole moments, 
they also proved that they completely characterize the 
isolated horizon geometry, that a unique isolated horizon 
(up to diffcomorphism) can be constructed from given 
multipole moments. 

A few years later, Schnetter, Krishnan, and Beyer [3l[ 
were the first to apply this multipole moment formalism 
in numerically generated dynamical spacetimes. Their 
work was intended as a wide overview of the use of the 
dynamical horizon formalism in interpreting numerical 
relativity simulations; for them, multipole moments were 
just one of many points of discussion. They applied the 
formalism of [301 ] in an essentially unmodified form. Be- 
cause this construction is restricted to axisymmetric hori- 
zons, the authors of (3l| focused attention on an axisym- 
metric black hole merger. 

Another application of this method appeared in 
Ref. [32j], a paper presenting methods to solve for con- 
formally curved initial data sets. As one might expect, 
when solving for a fully stationary single-black-hole ini- 
tial data set, the result is a slice of the Kerr spacetime, 
a fact that the authors confirm using the multipole con- 
struction of Ref. [3(J ■ 

Quite recently, another paper appeared [33| which in- 
troduced a novel scheme for computing multipole mo- 
ments indirectly, from surface integrals of various powers 
of the curvature. This new method is still restricted to 
axisymmetric horizons, but it avoids the need to explic- 
itly find the axisymmetry, and could markedly improve 
accuracy in cases where it can be used. 

Here, we take a slightly different approach. Rather 
than directly applying the methods of Ref. [30l ] in an ax- 
isymmetric merger, we modify the method, in a manner 
briefly suggested by its authors, so that it can be applied 
without the requirement of axisymmetry. Whereas the 
original method in Ref. [30j involved a preferred coordi- 
nate system on the axisymmetric horizon, in which spher- 
ical harmonic projections could be taken, we choose to 
project the relevant quantities against spectrally-defined 
spherical harmonics. Such harmonics are invariantly de- 
fined on any given topological sphere endowed with in- 
trinsic geometry, as eigenfunctions of geometric opera- 
tors, such as the one mentioned above relevant to the 
computation of spin angular momentum. Extra struc- 
ture, such as the axisymmetry that provides the preferred 
coordinate system of Ref. [30(, is not necessary. While 
the continuum eigenvalue problems that define these har- 
monics would complicate analytical treatments, they are 
quite straightforward to solve numerically. 



In section[TT]we introduce the details of this method, in 
particular the eigenvalue problems used to define spheri- 
cal harmonics on deformed spheres. In section iHll we in- 
vestigate one of the simplest applications of current phys- 
ical relevance. This is the question of the final remnant 
of a numerical merger of two vacuum black holes. While 
the general expectation is that the remnant of such merg- 
ers will generically be a Kerr black hole, relatively little 
effort has gone into a detailed investigation of whether 
this is actually the case. This question is of relevance to 
the status of black hole uniqueness, whose rigorous proof 
still involves certain analyticity assumptions [3~H ]. It is 
also related to the question of stability of the Kerr so- 
lution, which has so far been proven only for individual 
modes of linear perturbations (35j . Even if we fully ac- 
cept the expectation that general relativity must force the 
remnant of a black hole merger to be Kerr, the detailed 
recovery of the Kerr solution at late times, in as gauge- 
invariant a manner as possible, provides at the very least 
a stringent and physically-relevant code test. In Ref. [3r| . 
Campanelli ct al. demonstrated that a particular black 
hole merger simulation approaches Pctrov type D in a 
certain sense at late times, and carries no NUT charge. 
This fact largely confirms that their merger produces a 
Kerr geometry. One advantage of their approach is that 
it is fully local, that one can investigate the approach to 
Kerr geometry throughout the spatial slices, rather than 
simply on the horizon as we do here. In a followup to the 
current paper, we intend to repeat many of the methods 
of Ref. [36] on the datasets discussed in Sec. Mil Here we 
focus on multipole moments partly as a complementary 
method of black hole characterization, but also because 
these moments are of interest in their own right, as tools 
for probing the physics of tidal structure in strong-field 
gravity. 

II. GENERALIZED SPHERICAL HARMONICS 

The definitions given in Ref. [3(| for the mass and 
current multipoles on isolated horizons are very simple 
spherical harmonic projections of quantities related to 
the intrinsic and extrinsic geometry of the apparent hori- 
zon 2 in spacetime. 

I a := j> y a R dA, (2) 

L a := j>y*u B dA. (3) 

Here, dA refers to the metric volume element on the ap- 
parent horizon, R is its intrinsic scalar curvature (not to 



2 In the case of isolated horizons, the surfaces of interest are ar- 
bitrary two-dimensional spacelike slices of the three-dimensional 
null isolated horizon. In the case of dynamical horizons, the 
two-surfaces of interest are the apparent horizons into which the 
dynamical horizon is naturally foliated. 
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be confused with the Ricci scalar of the full spacetime, or 
of the spatial slice, or of the horizon worldtube), and lua 
is a connection on the normal bundle of the two-surface, 
which is conveniently written in terms of the two future- 
directed null normals, I and n: 

:= e A n v Wv/, (4) 

where ( 4 )\7 is the metric-compatible torsion- free space- 
time covariant derivative, and {e^} are basis vectors tan- 
gent to the two-surface. Throughout this paper, capital 
latin letters will index this two-dimensional tangent bun- 
dle. The null normals I and n are, as usual, normalized 
such that l-n = — 1. In most numerical papers and codes, 
u>a is written and computed in terms of the extrinsic cur- 
vature of the spatial slice. Here, we will refer to the I a as 
the mass multipoles and the L a as the current multipoles, 
though as noted in Ref. [3(| extra factors involving hori- 
zon areas and quasilocal spins must be included if one 
wishes to make them dimensionally consistent with the 
standard definitions of these quantities. 

The objects y a appearing in © and y A appearing in 
([3]) are scalar and vector spherical harmonics, respec- 
tively. It is in the definition of these harmonics that 
the breaking of axisymmetry has the most immediately- 
apparent cost, and therefore where the work in this pa- 
per will depart most strongly from the construction in 
Ref. [3(| ■ In that paper, attention is focused on the case 
of axisymmetric isolated horizons. Axisymmetry pro- 
vides a natural coordinate system on the apparent hori- 
zon, so the spherical harmonics used in Ref. [30| are the 
standard (m = 0) ones of spherical coordinates, applied 
in this canonical coordinate system. In other words, they 
are eigenfunctions not of the geometric Laplacian on the 
apparent horizon, but rather of the Laplacian of a metric 
sphere in these coordinates. There is nothing inherently 
wrong with such a choice in axisymmetry, in fact it pro- 
vides certain benefits in that context 3 , but for purposes of 
strongly dynamical, strongly non-axisymmetric systems 
a more general approach is called for. 

A. Scalar Spherical Harmonics 

Our approach will be to define the spherical harmonics 
spectrally, as eigenfunctions of the geometric Laplacian 
operator (or certain generalizations thereof) on the ap- 
parent horizon surface. In other words, our scalar spheri- 
cal harmonics are taken to be the functions y a that satisfy 
the equation 

&y a = \ a )Va (5) 



3 One such benefit is that the mass dipole moment always turns 
out to be zero. In other words, their construction guarantees 
that one is in a "center of mass frame." In general, this may 
not hold in our construction, though we have not yet seen an 
example where it fails. 



for some constant, X/ a y The function y a is defined only 
on the apparent horizon, and A is the intrinsic Laplacian 
of the apparent horizon, A := g AB V a^ b- The letter a 
is a label for the various solutions to the eigenproblem. 

Because the Laplacian in ([5]) reduces to the standard 
spherical Laplacian when the surface becomes a metric 
sphere, the functions y a reduce to the standard spherical 
harmonics in that special case as well. However, this is 
not the only self-adjoint operator with this property. For 
example, we can consider the problem: 

Ay a + qRy a = \ a )Va, (6) 

where R is again the intrinsic scalar curvature of the sur- 
face and q is a numerical parameter. In the case of a 
metric sphere, where R is constant, the second term on 
the left side does not alter the eigenfunctions, it merely 
increases each eigenvalue. This eigenproblem, therefore, 
can again be considered to define a reasonable general- 
ization of coordinate spherical harmonics. However, on a 
deformed sphere, where R is not constant, these general- 
ized spherical harmonics will no longer agree with those 
defined by ([5]). To fix this arbitrariness, and since we see 
no particular reason to prefer any other value for q, we 
choose q = 0, in other words the problem in Eq. ([5]), to 
define our scalar spherical harmonics. In the case of vec- 
tor spherical harmonics, we will see a geometrical reason 
to prefer a particular value for an analogous parameter. 



B. Vector Spherical Harmonics 

We will take our generalized vector spherical harmonics 
to be tangent to the surface, in which case they can be 
written in terms of gradients of two scalar potentials: 

y A = V A w a +e AB V B z a . (7) 

Here V is the torsion-free metric-compatible derivative 
on the apparent horizon, and cab is the Levi-Civita ten- 
sor on it. To consider the importance of these two po- 
tentials we should investigate the one-form lja against 
which the vector spherical harmonics will be projected. 
In Eq. (|4]), the future-directed null vectors I and n are 
orthogonal to the apparent horizon and normalized rel- 
ative to one another by the standard Newman-Penrose 
condition n ■ £ = — 1, but are otherwise free. One can 
arbitrarily scale the I vector at the cost of inversely scal- 
ing the n vector. This "boost freedom" is a standard 
gauge degree of freedom in the dynamical horizon for- 
malism. The dynamical horizon worldtube carries with 
it a preferred slicing into apparent horizons, but this slic- 
ing is only of the dynamical horizon itself. There is no 
preferred way of extending this slicing into the ambient 
spacetime. If we wish for our horizon multipoles to be 
independent of this gauge freedom, then we must choose 
harmonics that project out only the gauge invariant part 

of LUA- 
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From Eq. (HJ, it is apparent that a boost, £ n- al, 
fi M> a~ 1 ri, will add a pure gradient to uja'- 

uj a h-> w A - Va log(a). (8) 

Any part of that is a pure gradient is therefore en- 
tirely due to boost gauge, in the sense that it can be 
transformed away by an appropriate boost. Vector spher- 
ical harmonics of the form y A = S7 A w a will pick up this 
gauge-dependent information in the integral ([3]) , however 
vector spherical harmonics of the form y A = e AB V ' sz a 
will not. We therefore restrict all attention to vector 
spherical harmonics of this latter form. 

We now need a rule to define the potential functions 
z a that appear in these vector spherical harmonics. In 
the case of a metric sphere, the obvious choice is that 
they be the scalar spherical harmonics. As in the pre- 
vious subsection, there are many ways to generalize the 
spherical harmonics of the metric sphere. For the current 
purposes, there is reason to prefer a somewhat compli- 
cated fourth-order generalized eigenproblem: 

A 2 z a + V A (R\7 A z a ) = \ {a) Az a . (9) 

This generalized eigenproblem also defines the potentials 
for the approximate Killing vector fields used for comput- 
ing spin angular momentum in Ref. [22[ • For this reason, 
when this problem is used to define the vector spherical 
harmonics, the current dipole moment of the horizon is 
identical to the quasilocal spin defined there, a quantity 
that itself reduces on axisymmetric isolated horizons to 
the quasilocal spin defined by Hamiltonian methods [24| . 
In Ref. [30j], the agreement of the current dipole with 
the spin is cited as a reason to prefer using coordinate 
harmonics in a canonical coordinate system rather than 
spectrally-defined harmonics. There it was assumed that 
such harmonics would be simple eigenfunctions of the 
Laplacian, like the scalar spherical harmonics of the pre- 
vious subsection, in which case the current dipole would 
not agree with the standard spin angular momentum. We 
have averted this situation simply by choosing a better 
operator. 

We should also note that when vector spherical har- 
monics are chosen in this way, we are assured that there 
will be no current monopolc moment. This fact can be 
viewed in a number of related ways. On the simplest 
level, there is the fact that when the vector spherical 
harmonics are defined to be of the form y A = e AS V 'sz a , 
then a potential of the form z a = const, cannot define a 
(normalizable) vector spherical harmonic. In some sense, 
z a = const, can be viewed as a solution to Eq. §§§ with ar- 
bitrary eigenvalue, but it is not a well-behaved solution. 
The generalized eigenproblem is technically singular in 
function spaces that include constants [13] , meaning that 
well-behaved solutions cannot be found unless the func- 
tion space is restricted to, for example, functions with 
zero average over the sphere, a condition which removes 
all nonzero constants from consideration. 

Another way of looking at this, which helps to eluci- 
date the relationship between the mass and current mul- 



tipoles, is that when the vector spherical harmonics are 
defined in this way, an integration by parts allows the 
current moments to be written as: 

L a := <j> z a n dA, (10) 

where il := € AB \7a<^b can be interpreted geometrically 
as a scalar curvature of the normal bundle of the two- 
dimensional surface in four-dimensional spacetime. The 
current moments thus represent for the extrinsic geom- 
etry of the apparent horizon what the mass moments 
represent for its intrinsic geometry. 

The complex combination of these two curvatures, R + 
iCl, is sometimes called the complex curvature of the two- 
surface embedding. As is briefly described in Sec. 4.14 
of Ref. [HI , the vanishing of the current monopole mo- 
ment can be understood geometrically in this context as 
a result of the generalization of the Gauss-Bonnet theo- 
rem to the lorentzian normal bundle. The integral of any 
constant multiple of is a topological invariant, just like 
that of R, but because the gauge group on the normal 
bundle is topologically trivial, this invariant must always 
vanish. 

One final point to note, with regard to both the scalar 
and the vector harmonics, is that of normalization. So- 
lutions of the cigenproblems in Eqs. (J5)) and (|9]) are de- 
termined only up to constant multiplicative factors 4 . We 
fix these factors with an integral normalization condition. 
The condition imposed on scalar spherical harmonics is: 

I (y a f dA = 1. (11) 

On metric spheres in Euclidean space, this reduces to 
the standard normalization condition for scalar spherical 
harmonics (up to a factor of areal radius). 

For vector spherical harmonics, the normalization con- 
dition wc impose is: 

jg A BvivldA = \. (12) 

This differs slightly from the standard normalization con- 
dition for axial vector spherical harmonics in euclidean 
space, which involves an extra factor of £{£ + 1), but be- 
cause the generalization of the index £ is not an integer, 
but rather a function of the eigenvalue A( Q ), we simply 
leave this factor out. 

Normalization conditions like those above still don't 
determine a sign for the spherical harmonics. This sign 
ambiguity translates directly into a sign ambiguity for 
the multipoles. We fix the sign with the condition that 
the values of the multipolc moments be nonnegative. 



4 In fact, solutions of Eq. (0 arc determined only up to constant 
multiplicative and additive factors, however additive constants 
have no effect on the multipoles due to the vanishing of the cur- 
rent monopole. 
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In summary, the approach we take to defining mul- 
tipoles in numerical simulations begins with finding so- 
lutions y a and z a of the eigcnproblcms ([5]) and ([9]) de- 
fined on the apparent horizon two-surface. From the po- 
tential z a , vector spherical harmonics are computed as 
Va : ~ £ AB ^ BZ a - The harmonics are normalized with 
the conditions in Eqns. (fTTj) and (fT2]l. Then, the surface 
integrals in Eqns. @ and © arc computed to be the 
multipoles. 



III. NUMERICAL RESULTS 

The immediate purpose of this mathematical machin- 
ery is to investigate the remnant of a binary black hole 
merger. There is a very large space of physically-relevant 
mergers worth investigating, including variations in the 
initial mass ratio, eccentricity, spin magnitudes, and spin 
directions. For this paper we will focus on a very simple 
case: the merger of a non-eccentric binary of equal mass, 
nonspinning black holes. This data set is discussed in 
detail in Ref. [3] , which briefly notes the fact that two in- 
dependent measures of the final spin agree to well within 
their expected numerical errors. This claim can be con- 
sidered a first indication that the tidal structure of the 
quiescent black hole is that of Kerr, as this is the case in 
which these two measures of spin are designed to agree. 
Our goal now is to present the rest of the tidal informa- 
tion, to the extent that it can be resolved in the code, 
to strengthen the case that the final remnant is a Kerr 
black hole. 

The code used to compute these multipoles is a part of 
the Spectral Einstein Code (SpEC) developed and main- 
tained by the Caltech and Cornell Numerical Relativity 
groups, particularly Lawrence E. Kidder, Harald P. Pfcif- 
fer, and Mark A. Scheel. Once an apparent horizon has 
been found, using the method described in Ref. [39| , the 
code interpolates all relevant data to a pseudospcctral 
grid on that surface. Because this grid is pseudospectral, 
the code can automatically transform any smooth func- 
tion on the apparent horizon into a truncated expansion 
in coordinate spherical harmonics. This expansion, in- 
serted into Eq. ([5]) or ((9]), provides a finite-dimensional 
matrix eigenproblem (or generalized eigenproblem in the 
latter case) which is solved using the LAPACK routine 
dggev. We emphasize that while this construction in- 
volves coordinate spherical harmonics, they are only used 
to supply the numerical discretization, so they disap- 
pear in the continuum limit. The geometrical spheri- 
cal harmonics y a and that define the multipoles are, 
apart from numerical truncation error, uniquely defined 
on any given two-sphere (up to possible degeneracy in 
the eigenspaces). 

The information that we can assess includes not only 
the values of the multipole moments defined in Eqns. ([2]) 
and ©, but also the spectrum of eigenvalues in Eqns. ([5]) 
and A particular motivation for investigating eigen- 
values of geometric operators is that they provide an in- 



dication of symmetries in the horizon. As is familiar from 
elementary quantum mechanics, a symmetry in an opera- 
tor leads to degeneracies in its eigenspaces. The converse 
is not necessarily true, but on an intuitive level we may 
interpret degeneracies in the eigenspectrum as indicators 
of possible symmetry. 

This is an interesting tool for the study of this particu- 
lar problem, because in the ringdown after a nonspinning 
black hole merger there is a transition from one axis of 
symmetry to another. Immediately after the formation 
of a common apparent horizon, one intuitively expects 
this horizon to be "peanut shaped," with an axis of ap- 
proximate symmetry 5 along a line connecting the pre- 
vious two individual apparent horizons. After the ring- 
down is complete, one would expect a single black hole 
with symmetry about the axis of the initial orbital an- 
gular momentum. This breaking and forming of sym- 
metries is demonstrated in Fig. [1] The figure presents 
the three eigenvalues of the horizon Laplacian associated 
with harmonics that would settle to the £ — 1 spheri- 
cal harmonics if the horizon were to become metrically 
spherical. Two of these curves overlap at early times, 
a degeneracy due to the approximate axisymmetry of 
the initial "peanut" shape. As this symmetry is bro- 
ken during the ringdown, the degeneracy breaks and one 
eigenvalue eventually joins up with the third eigenvalue, 
demonstrating the eventual axisymmetry about the spin 
direction. 6 

With the next five eigenvalues, in Fig. [2J we see the 
pattern again. Again, modes are nearly degenerate at 
early times, but split off during the ringdown and re- 
connect as the quiescent symmetry is approached. Note, 
however, that one degeneracy at early times is quite vis- 
ibly broken. This may be due, on an intuitive level, to 
the tidal interaction of the two black holes during mspi- 
ral, with shifted phase due to horizon viscosity [40l442l | . 
Because such tidal interaction is a quadrupolar effect, it 
would make sense for it to be less visible in the dipolar 
information of Fig. [T] 

Degeneracies in the other eigenproblem, Eq. ©, give 
a similar picture of the breaking and reforming of sym- 
metries, but this problem gives an even more compelling 
picture of the relationship between symmetries and de- 



In the case studied here, this axisymmetry would only be ap- 
proximate, as tidal bulges would be expected to phase-shift due 
to horizon viscosity during the inspiral. In the case of a direct 
head-on collision of nonspinning holes, this axisymmetry would 
be exact, and would be preserved even through the ringdown. 
All figures in this paper give quantities computed in code units 
evaluated with respect to coordinate time which is also expressed 
in code units. For context, the final black hole described in 
Figs. [l]{9] has horizon mass M pa 1.98 in these code units, where 
this mass is defined by the Christodoulou formula M 2 = Mf rr + 
J 2 /(4Af a 2 rr ) where Mi rr = A/(16n) is the irreducible mass and 
J is the quasilocal spin angular momentum defined in Appendix 
A of Ref. [52] • In the simulation presented in Figs. IT0U13I the 
value of this final mass in code units is M « 2.56. 
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FIG. 1: Absolute values of the lowest three (dipole) non- 
trivial eigenvalues of the Laplacian on the dynamical horizon 
during the ringdown. The breaking away of the red (dotted) 
curve from the blue (solid) curve at early times is due to the 
breaking of the initial (approximate) "peanut" axisymmetry. 
The joining of this curve onto the black (dashed) curve at late 
times is due to the late-term axisymmetry of the final Kerr 
horizon. The horizontal gray lines represesent the expected 
eigenvalues on a Kerr horizon with mass and spin equal to the 
final measured values in the simulation. Thus the convergence 
of the eigenvalues to these lines demonstrates the approach to 
Kerr geometry. 
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FIG. 2: Absolute values of the next five (quadrupole) eigen- 
values of the Laplacian on the dynamical horizon during the 
ringdown. As in Fig.[T] the curves are paired up at early times, 
split their degeneracies, and connect in a different pairing at 
late times, again indicating transition from one axis of sym- 
metry to another. The fact that one of these degenerate pairs 
at early times is visibly nondegenerate indicates imperfection 
in the "peanut" axisymmetry intuitively due to phase offset 
tidal bulges built up during the inspiral. 



FIG. 3: Absolute values of the lowest three (dipole) nontrivial 
eigenvalues of the generalized eigenvalue problem in Eq. (J9j) 
on the dynamical horizon during the ringdown. The vertical 
axis is now scaled logarithmically to better show the approach 
of the smallest eigenvalue to zero. As argued in Appendix A 
of Ref. [22|. the vanishing of this smallest eigenvalue is di- 
rect evidence of a rotational symmetry of the intrinsic surface 
geometry, so this figure provides a clear picture not only of 
the symmetry transition itself, but also of the relationship be- 
tween symmetries and degeneracies. In particular, the cross- 
ing of the red and black (dotted and solid) curves can be seen 
as an example of "accidental" degeneracy, degeneracy that is 
not necessitated by a symmetry of the operators. 

generacies. The original motivation of Eq. ©, as de- 
scribed in Appendix A of Ref. (22j . was to construct, in 
a sense, the closest possible approximation of an axial 
symmetry on a horizon that may not have any true sym- 
metries at all. One can easily show that the value of 
a given eigenvalue is proportional to the integral of the 
square of the residual in Killing's equation for the asso- 
ciated "approximate Killing vector" field. Thus, when 
there is a true symmetry, and therefore a true Killing 
vector field, one of the eigenvalues of this problem will 
equal zero. So from plots of the eigenvalues of Eq. © 
we can see the breaking and forming of symmetries both 
indirectly, through degeneracies of the eigenspaces, and 
directly, through the value of the lowest eigenvalue. Fig- 
ure [3] shows the three lowest eigenvalues of this problem. 
As noted at the end of Sec. Ill B[ there are no monopole 
harmonics at all for this problem, so these are the three 
harmonics that would reduce to the I = 1 harmonics if 
the horizon approached a metric sphere. The vertical 
axis of the figure is now logarithmically scaled, to show 
the approach of the smallest eigenvalue to zero both at 
early and late times. 

These figures also provide a quantitative picture of the 
intrinsic geometry of the apparent horizon, and its ap- 
proach at late times to the geometry of a slice of the 
Kerr horizon. The horizontal lines in Figs. [1] [3] represent 
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FIG. 4: Relaxation of the three excited quadrupole moments 
to expected Kerr values. Two of the five possible quadrupole 
moments must vanish due to the reflection symmetry in the 
problem, and indeed, their computed values are small enough 
to be considered zero to within ordinary numerical errors. Of 
the remaining three multipoles, two fall exponentially toward 
the level of numerical truncation, and the third quickly settles 
to the expected value for a Kerr black hole of the same final 
mass and spin. This expected value is shown in the thick 
horizontal gray line, which for most of the simulation over- 
laps the dashed blue curve. In this and later figures, curves 
are labeled by the "effective" I values of their corresponding 
spherical harmonic at late times. That is, £ e s is defined from 
the value of the eigenvalue A at late times by the equation 
A = — l c s(£c« + 1)A' 2 , where r is the areal radius of the hori- 
zon, again measured after the hole settles down. 



the expected values for these eigenvalues on a Kerr hori- 
zon of the same mass and spin as is measured at very late 
times in the simulation. Note that this spin is guaranteed 
to be identical to the late time current dipole moment on 
the horizon (this is the main motivation for Eq. ((9|)), so 
agreement of the current dipole with the "expected Kerr 
value" is trivial, however the consistency of all other mul- 
tipoles, as well as these eigenspectra, present a nontrivial 
demonstration that the quiescent hole is Kerr. 

Figures 0H5] present the behavior of the multipole mo- 
ments. In Fig. 31 the three excited quadrupole moments 
are shown (the other two vanish as demanded by reflec- 
tion symmetry). One moment starts out relatively small 
and grows to take the value expected for a Kerr black 
hole. The other two fall exponentially toward zero, un- 
til reaching the level of numerical truncation. Figure [5] 
shows the convergence of this floor of numerical error for 
three values of the resolution of the numerical simula- 
tion. On all three simulations, the horizon finder and 
eigenvalue solver are run at the maximum relevant reso- 
lution, essentially the same as the angular resolution of 
the original simulations. Figures [MS are analogous to 
Fig. 21 showing higher-order multipole moments. Again, 
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FIG. 5: A particular quadrupole moment from Fig. U shown 
for three different resolutions. The order of the pseudospectral 
angular discretization is given by L = 18, 20, 22 respectively, 
L representing the maximum £-value of coordinate spherical 
harmonics used to discretize the problem. At late times, the 
exponential falloff halts, but the level where this occurs con- 
verges exponentially toward zero as L is increased. These 
nonzero values can therefore be attributed to standard trun- 
cation error. 



all multipoles allowed by the reflection symmetry of the 
problem are excited near the moment of merger, but in 
each case a single moment rises to the expected value 
for a Kerr black hole of the measured final mass and 
spin, while all other multipoles decay exponentially to- 
ward zero before stopping due to numerical truncation. 

The fact that those multipoles that decay to zero do 
so exponentially raises the question of whether this decay 
can be attributed to quasinormal ringing. The answer to 
this question is clouded by a few subtleties. For one, 
while the multipoles do appear to oscillate within an ex- 
ponential envelope, this oscillation does not appear to 
be even approximately periodic, and at any rate occurs 
on a much longer timescalc (compared to the exponen- 
tial decay timescale) than the oscillations associated with 
quasinormal ringing. There is an intuitive explanation 
for this. Because the multipoles arc defined with respect 
to spherical harmonics that are fixed by the intrinsic ge- 
ometry of the horizon, changes in horizon geometry will 
cause changes in these harmonics. In particular, if the 
major part of the perturbation from Kerr geometry is 
a nonaxisymmetric bulge that rotates around the spin 
axis, then the spherical harmonics will be dragged along 
with this bulge. Intuitively, an ideal "£ = 2, m = 2" 
bulge would be expected to drag the spherical harmonics 
into corotation with it, so the multipole representing this 
bulge would be expected to fall off as a pure exponen- 
tial, with no oscillation. In reality, the situation is more 
complicated, presumably due in part to the existence of 
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FIG. 6: Relaxation of the five excited hexadecupole mass 
moments to expected Kerr values. Four of the nine possible 
moments vanish due to the reflection symmetry in the prob- 
lem. Of the remaining five, four fall to zero exponentially, and 
the other rises to the expected value for a Kerr black hole. 
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FIG. 7: Of the seven current octupole moments, only three 
are allowed by the reflection symmetry to be excited. Again, 
two fall exponentially toward the level of numerical truncation 
error, and the other exponentially approaches its expected 
value for a Kerr hole. The quantity £ e g labeling the curves 
is again computed from the eigenvalue at late times, but the 
definition is now given by the equation A = — £ c b(£ c h + l)/r 2 + 
2/r 2 , taken from the special case of Eq. (J5J) on a round sphere. 



higher multipolar structure, and in part due to the even- 
tual approach to axisymmetry, causing degeneracies in 
the eigcnproblems to be broken at the numerical level 
rather than at the analytical one. 

Properly "unwinding" this rotation of the harmonics 
would amount to a partial fixing of angular coordinates. 



FIG. 8: Of the eleven current 32-pole moments, only five are 
allowed by the reflection symmetry to be excited. Again, four 
fall exponentially toward the level of numerical truncation 
error, and the other exponentially approaches its expected 
value for a Kerr hole. The quantity £ c h labeling the curves is 
defined as in Fig. [7] 



There may be sensible ways to do this, but we consider 
this somewhat oustide the scope of the current research, 
so instead we choose to ignore the oscillatory behavior, by 
focusing on the quadratic sums of multipoles associated 
with nearly-degenerate eigenspaces. In particular, the 
two exponentially-decaying curves in Fig. 0] are associ- 
ated with such an asymptotically degenerate eigenspace, 
and can intuitively be interpreted as real and imaginary 
parts of the u £ = 2, m = 2" multipole. Their quadratic 
sum can therefore be viewed as the overall "magnitude" 
of the quadrupolar part of this rotating bulge, and would 
be expected to fall off exponentially in time without os- 
cillation. Figure [§] shows the value of this quadratic sum 
as well as the falloff rate expected from perturbation the- 
ory. We used the method due to Leaver [43| to compute 
quasinormal frequencies in terms of the roots of two cou- 
pled complex continued fractions. For the mass and spin, 
we used values computed at late times on the horizon and 
reported in Ref. 

The remarkably fine agreement between this quadratic 
sum of multipole moments and the expected exponential 
falloff of the dominant "£ — 2, m = 2" quasinormal mode 
in perturbations of the Kerr geometry makes it tempting 
to try to pick other quasinormal ringing modes out of 
the data. This however would be a somewhat nontrivial 
undertaking. For one thing, all multipoles defined in our 
formalism are computed from data directly on the hori- 
zon. Different radial modes of black hole perturbations 
would be directly superposed in the multipoles. Further 
complicating matters in the case of a Kerr hole, the an- 
gular dependence of the quasinormal modes are not given 
by pure spherical harmonics, as defined here, but as solu- 
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FIG. 9: The quadratic sum Q(t) := ^qf(t) + q%(t) of the two 
moments qi and qi that fall exponentially to zero in Fig. 0] 
The thick gray line shows the expected falloff of this quantity 
in perturbation theory. 



tions of the "angular equation" of the Teukolsky formal- 
ism (Eq. (4.10) of Ref. [4J])- Thus all multipoles would be 
expected to project out components of all radial and an- 
gular quasinormal modes, rather than cleanly projecting 
out one at a time. What is seen in Fig. 13 as exponential 
decay is actually just the dominant term of a multiexpo- 
ncntial expansion. The problem of fitting data to a sum 
of exponentials is famously ill posed, so any effort to pick 
out higher-order ringing modes from this data would be 
quite ambitious, if possible at all. 

Many of these complications disappear if the final black 
hole is nonspinning. In that case, the quasinormal modes 
should have the angular dependence of pure spherical 
harmonics, so multipoles of a given order can be ex- 
pected to project out modes of the same order (though, 
again, multiple radial modes would be expected to over- 
lap). Also, if there arc enough degrees of reflection sym- 
metry to forbid the rotation of the spherical harmonics 
described above, one might hope to recover not only the 
exponential falloff rates of different quasinormal modes, 
but also their frequencies of oscillation. 

Figures [TUHTSl demonstrate this recovery. The data 
used here are from the ringdown after the collision of 
two black holes of nonzero antialigncd spin (and there- 
fore zero total angular momentum) starting from rest. 
This is a simple test case that can be used for studying 
black hole kicks, and the particular simulation studied 
here will be presented in great detail for that goal in an 
upcoming paper [45[. For the present purposes, the im- 
portant points are that the final state is nonspinning, and 
that two orthogonal planes of reflection symmetry (the 
coordinate x = and z = planes, with the final kick 
being in the y direction) hold the harmonics in place, in 
that they remain symmetric or antisymmetric under the 



action of the reflection symmetries. Figures analogous 
to Figs. [TU1 - H51 for the current multipoles show similar 
agreement, but are omitted here because they look es- 
sentially the same. Incidentally, similar (though less de- 
tailed) agreement with quasinormal ringing frequencies 
was noted in the oscillation of the area of spatial slices 
of the event horizon in other recent simulations using the 
SpEC code [Hj|. 

One subtlety with Figs. [TU1 flUl must be noted. At late 
times, when these moments reach levels on the order of 
1CP 9 , the data become quite noisy. This is a result of nu- 
merical errors, particularly the truncation error of the an- 
gular discretization. However it appears that the numeri- 
cal data in Figs. [Tl~r[T3l continue to decrease as the simula- 
tion goes on. This is an artefact of the manner in which 
the data are extracted. The numerical code computes 
essentially as many multipoles as there are grid points 
on the interpolated apparent horizon. These multipoles 
must be ordered in some way. The most obvious ordering 
is provided by the eigenvalues of the spherical harmonics. 
However such an ordering is not effective when families of 
eigenspaces are nearly degenerate, as particularly in the 
case of ringdown to a Schwarzschild black hole. To pick 
out particular eigenvalues in this quasinormal ringdown 
phase, we employ a simple post-processing script that 
chooses, at each timestep, the particular multipole mo- 
ment that has value closest to a "prototype" value taken 
from the perturbation theory results shown in the red 
curves in these figures. Simultaneously, the script checks 
that any chosen multipole corresponds to an eigenvalue 
which lies within a certain range, so that the multipole is 
assured to have the proper value. After this search- 
ing is carried out, we check the eigenvalue, as a function 
of time, corresponding to the chosen multipole, to ensure 
that it is smooth and therefore that the procedure has 
chosen a consistent multipole and eigenvalue dataset 7 . 
This procedure nicely and unambigously recovers physi- 
cal perturbations during most of the ringdown. However 
at late times the perturbation is small enough that the 
ordering ambiguity is particularly strong. At late times, 
the script chooses the moment closest to the prototype 
value, out of the many that are oscillating quickly at 
small values due to numerical error, yet all eigenvalues in 
the given range are essentially the same, so the smooth- 
ness of the eigenvalue is no longer an effective tool to dis- 
tinguish the correct moment from the others of the same 
£. The matching of the numerical data to the prototype 
function is therefore given more weight than it deserves, 
and the data, though clearly flooded with numerical er- 



7 Immediately after the formation of the common horizon, 
when deviations from the expectations of linearized theory are 
strongest, this script can again have trouble finding consistent 
multipole and eigenvalue datascts. For this reason some data 
are omitted from the beginning of Figure [T2l as nonsmoothness 
of the eigenvalues showed that the chosen multipoles did not 
represent a consistent dataset at very early times. 
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FIG. 10: Ringing of the dominant mass quadrupole mo- 
ment as the product of a nonaxisymmetric head-on merger 
settles down to a Schwarzschild black hole. The dashed blue 
curve represents numerical data, and the solid red curve is 
an exponentially-damped sinusoid, with frequency and damp- 
ing corresponding to quadrupole gravitational quasinormal 
modes of a Schwarzschild black hole [4j|. An arbitrary con- 
stant amplitude scaling and phase have been applied to the 
red curve, by eye. 



ror, continue to fall off exponentially in time. Figure [TU] 
is an exception to this behavior. In that particular case, 
a method involving matching the spherical harmonics to 
coordinate spherical harmonics was able to unambigu- 
ously pick out the "correct" harmonic. For higher multi- 
poles that method failed, apparently due to the rotation 
ambiguity of the coordinate spherical harmonics them- 
selves. 

The quality of the agreement with standard quasi- 
normal ringing frequencies, both in the approach to 
Schwarzschild geometry in Figs. [TOHT51 and in the ap- 
proach to Kerr geometry in Fig. initially came as quite 
a surprise, considering that a major motivation for this 
project is a healthy skepticism for the quality of the coor- 
dinates in numerical simulations. While the slicings used 
here are not arbitrary — both simulations employ har- 
monic slicings during the ringdown, as do conventional 
treatments of black hole perturbation theory — they are 
nonetheless different harmonic slicings than those in con- 
ventional perturbation theory, because the ones used in 
our simulations arc horizon-penetrating. One might ask, 
then, how the numerical code knows to settle on a har- 
monic slicing in which these frequencies come out as ex- 
pected. The answer lies in the approach to stationarity. 
At late times the simulations develop approximate sta- 
tionary Killing vector fields £, and the coordinate compo- 
nents of the spacetime metric tensor asymptote to con- 
stant values, meaning that the coordinates adapt them- 
selves to the symmetry such that dt — > £. This turns the 
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FIG. 11: Ringing of the dominant mass octupole moment 
as the product of a nonaxisymmetric head-on merger settles 
down to a Schwarzschild black hole. 
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FIG. 12: Ringing of the dominant mass hexadecupole mo- 
ment as the product of a nonaxisymmetric head-on merger 
settles down to a Schwarzschild black hole. 



definition of the ringing frequencies into a geometrical 
statement: rather than saying df$ = —u) 2 $, one can say 
£ (£ (®)) = — w 3>. In other words, the frequencies come 
out right because the coordinates adapt themselves to the 
late-term stationarity. The process by which this adap- 
tation occurs is, to our knowledge, still an open question. 

At any rate, we must caution that this recovery of 
standard frequencies at late times should by no means 
be taken as license to overlook gauge ambiguity in nu- 
merical simulations. For example, it is quite tempting to 
associate the slight disagreement with perturbative re- 
sults immediately after merger with nonlinear dynamics, 
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FIG. 13: Ringing of the dominant mass 32-pole moment 
as the product of a nonaxisymmetric head-on merger settles 
down to a Schwarzschild black hole. 



however this disagreement could just as likely be due to 
the coordinates having not yet adapted to the approxi- 
mate stationarity, or to stationarity simply not existing 
to a sufficient approximation. All of these effects (and 
perhaps others) will have an impact on the ringing imme- 
diately after merger, and a detailed investigation of the 
nonlinear extension of quasinormal ringing would require 
(at least partially) slicing-invariant comparisons beyond 
the scope of the current work. For example, one might 
treat the ringing of one multipole as a "clock" by which 
to measure the frequencies of the other multipolcs. 



IV. DISCUSSION 

We have presented a definition of quasilocal source 
multipoles on dynamical horizons, adapted from that in 
Ref. [30j in such a way that it can be applied to hori- 
zons without axisymmctry, while preserving the agree- 
ment of the current dipole moment with the spin angular 
momentum defined by Hamiltonian methods 2J]. More 
precisely, the vector spherical harmonics used to project 
out current multipoles are constructed in such a way that 
the dipole moment is identical to the spin angular mo- 
mentum used in Ref. [22| . The key to this generalization 
is the definition of spherical harmonics as solutions to 
certain eigenvalue problems on the apparent horizon. 

We have also applied this formalism to demonstrate 
that in a detailed and partially gauge-invariant sense, the 
binary black hole merger described in Ref. [1] indeed set- 
tles to a Kerr black hole, at least in the neighborhood of 
the horizon. There are limits to the gauge independence 
of this statement. The work here depends heavily on 
the formalism of dynamical horizons [25| , which are de- 



pendent on the slicing of spacetime (or, from a different 
viewpoint, are themselves invariantly defined yet carry 
unique foliations into apparent horizons that are compat- 
ible with the foliation of spacetime only in certain time 
slicings). Use of a unique and invariantly defined horizon 
such as the event horizon may be of interest (and is possi- 
ble in the SpEC code [HI), however it would not alleviate 
the problem of slicing dependence, as a slicing must be 
chosen at some point to break the three-dimensional hori- 
zon worldtube into two-dimensional surfaces on which 
the spherical harmonic projections are taken. 

A demonstration along the same lines as discussed here 
has been carried out before j3l| . however the numerical 
results here are somewhat stronger, and our generaliza- 
tion of the formalism has allowed the consideration of a 
non-axisymmetric merger. 

Looking in detail at the ringdown of the multipoles, we 
have also recovered known quasinormal ringing frequen- 
cies. The dominant exponential damping timescale is 
recovered in the ringdown to Kerr geometry, and agrees 
with results from perturbation theory. Much more de- 
tailed results are found in the ringdown after a head-on 
collision leading to a Schwarzschild geometry, in which 
oscillation frequencies and damping timescales can be 
picked out mode by mode. 

In future work we intend to study the ringdown of these 
datasets (and possibly others) on a local level, using a 
variant of the method presented in Ref. [36| . 

As for the multipole moments themselves, various av- 
enues of investigation are open. The methods used here 
could be applied to study the tidal interaction of black 
holes during fully nonlinear binary inspiral and merger, 
inclu ding a full nonlinear generalization of certain re- 
sults [4 ll [42l] of black hole perturbation theory. As men- 
tioned in Ref. [30], quasilocal source multipoles might 
also be applicable in trying to find a generalization, to ex- 
act general relativity, of Einstein's celebrated quadrupolc 
formula. Related to this, one might hope to recover 
force laws at the quasilocal level, relating black hole kicks 
to products of multipoles, as is done in the asymptotic 
regime in Ref. (l7| . However, such an investigation would 
presumably require a satisfactory quasilocal definition of 
black hole linear momentum, which (if possible at all) 
appears to be beyond the realm of current understand- 
ing. 
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